Effect of CYP3A4 Inhibitors and Inducers on Pharmacokinetics and Pharmacodynamics of Saxagliptin and Active Metabolite M2 in Humans Using Physiological-Based Pharmacokinetic Combined DPP-4 Occupancy

We aimed to develop a physiological-based pharmacokinetic and dipepidyl peptidase 4 (DPP-4) occupancy model (PBPK-DO) characterized by two simultaneous simulations to predict pharmacokinetic (PK) and pharmacodynamic changes of saxagliptin and metabolite M2 in humans when coadministered with CYP3A4 inhibitors or inducers. Ketoconazole, delavirdine, and rifampicin were selected as a CYP3A4 competitive inhibitor, a time-dependent inhibitor, and an inducer, respectively. Here, we have successfully simulated PK profiles and DPP-4 occupancy profiles of saxagliptin in humans using the PBPK-DO model. Additionally, under the circumstance of actually measured values, predicted results were good and in line with observations, and all fold errors were below 2. The prediction results demonstrated that the oral dose of saxagliptin should be reduced to 2.5 mg when coadministrated with ketoconazole. The predictions also showed that although PK profiles of saxagliptin showed significant changes with delavirdine (AUC 1.5-fold increase) or rifampicin (AUC: a decrease to 0.19-fold) compared to those without inhibitors or inducers, occupancies of DPP-4 by saxagliptin were nearly unchanged, that is, the administration dose of saxagliptin need not adjust when there is coadministration with delavirdine or rifampicin.


INTRODUCTION
Patients undergoing multiple comorbidities are usually treated by complicated polypharmacy schemes and long-term administration, which means that they will be at the risk of drug-drug interactions (DDIs). In DDIs, one drug may affect the PK or PD behavior of another drug, thereby resulting in a series of side effects and even a withdrawal of approved pharmaceuticals from the market, for instance, mibefradil and terfenadine (Valicherla et al., 2019;Pelkonen et al., 2020). Therefore, prospective assessment of a potential risk of DDIs is significant within the pharmaceutical industry, such as guaranteeing the safety and reducing unnecessary consumption of new drugs. Cytochrome P450, which is the most comprehensive metabolizing enzyme in the gut and the liver, plays an important role in pharmacokinetic (PK) interaction-related DDIs through the interfering metabolism of most drugs (Dmitriev et al., 2019). Currently, PK-related DDIs mediated by CYP have become a research focus over the past decades (Conner et al., 2018;Peng et al., 2021). However, only focusing on PK-related DDIs is far from enough, and pharmacodynamic (PD)-related DDIs should receive more attention, particularly on dose adjustment for patients. Although PD-related DDIs recorded are roughly 1.9fold higher than PK-related DDIs (Spanakis et al., 2019), unfortunately, simulation study on PD-related DDIs is far less compared to those on PK-related DDIs.
Ketoconazole and delavirdine are strongly competitive and timedependent inhibitors (TDIs) against the CYP3A4 enzyme, and rifampicin is a strong CYP3A4 inducer. It is well known that ketoconazole and rifampicin are often recommended to assess potential DDIs for other drugs mainly metabolized by the CYP3A4 enzyme (Rytkönen et al., 2020;Xu et al., 2021). To our knowledge, owing to irreversible loss of P450 enzyme functions, TDIs are regarded to have a longer persistent time on the drug metabolic enzyme and consequently to cause clinically more significant DDIs in contrast to competitive inhibitors (Kosaka et al., 2017;Eng et al., 2020;Tanna et al., 2021). Therefore, ketoconazole, delavirdine, and rifampicin were chosen for potential DDI assessments.
Saxagliptin is an orally potent competitive DPP-4 inhibitor and utilized for the treatment of type 2 diabetes at a daily dose of 5 mg. 1 . Saxagliptin is metabolized extensively in humans via CYP3A4 to be converted into many metabolites, of which 5-hydroxy saxagliptin (M2) is the most major active metabolite. M2 in vivo is roughly in a 2-fold higher amount (44.1 versus 24.0%) and binds to DPP-4 with a ∼2-fold lower affinity than saxagliptin (1.3 versus 2.6 nm) 1 (Su et al., 2012). Hence, when coadministered with inhibitors or inducers of CYP3A4, significant effects of them on PK and PD of saxagliptin should be taken into account.
We aimed to develop a mathematical model to assess the dynamic effect of ketoconazole, delavirdine, and rifampicin on the PK and PD of saxagliptin tablets in humans, when coadministered. More precisely, the physiological-based pharmacokinetic (PBPK) model and the DPP-4 occupancy (DO) model were incorporated into a new mathematical model, termed as the PBPK-DO model, which enables the changes of PK and PD of saxagliptin and M2 in humans to be quantified simultaneously with coadministration of CYP3A4 inhibitors or inducers.

Data Collection
Clinical PK studies of ketoconazole, delavirdine, and rifampicin were taken from the published literature 2, 3 (Hanke et al., 2018), which could be used in their respective PBPK model establishment and verification. Clinical PK and PD studies of saxagliptin and metabolite M2 were collected from published data (Upreti et al., 2011), which could be used in the PBPK-DO model development and verification for saxagliptin and M2. Physicochemical properties of four drugs, binding kinetics of saxagliptin and M2, and physiological parameters in humans required in developing the PBPK-DO model were obtained from published scientific studies, and the built-in libraries of Gastroplus software were used in this study, including Berkeley Madonna (Version 10.2.8, Berkeley Madonna, Inc Albany, CA, United States) and ADMET Predictor (Version 9.0.0.0, Simulation Plus, Inc Lancaster, CA, United States).

Development of the PBPK-DO Model
A PBPK-DO model was developed to simulate PK and DPP-4 occupancy time profiles of saxagliptin and metabolite M2 simultaneously after oral coadministration with ketoconazole, delavirdine, and rifampicin. The PBPK-DO model was composed of three key simulation processes, that is, PK prediction of saxagliptin and M2, simulation of DPP-4 occupancy by saxagliptin and M2, and interaction prediction with inhibitors/ inducers of CYP3A4. In this PBPK-DO model, first, the simulation of saxagliptin and M2 concentration changes over time was enabled simultaneously by the PBPK model, which consisted of a stomach-gut compartment, a enterocytes compartment, a portal vein compartment, a blood compartment (arterial and venous blood), eliminating tissues (the liver, the kidney), non-eliminating tissues (adipose, the bone, the brain, the heart, muscle, the skin, and the spleen), and the lung. Next, the interactions between inhibitors/inducers and CYP3A4 were calculated through inhibition or inducing parameters (K i , K inact , E max , and EC 50 ), CYP3A4 expression amount, and free saxagliptin concentration in the gut and liver. Eventually, time courses of DPP-4 occupancy were characterized by two key rate constants of on-rate (k on ) and off-rate (k off ) combined with free drug concentration around the DPP-4. In addition, we assumed that DPP-4 was located in the venous blood compartment in the present model. The overall framework of the PBPK-DO model is represented in Figure 1.

Stomach-Gut Compartment
Assuming that the drug in the stomach was neither absorbed nor metabolized, the drug amount in the stomach (A 0 ) is only governed by the gastric emptying rate (K 0 ). The change in mass within the stomach is described as follows (Qian et al., 2019): The gut lumen consists of the duodenum, jejunum, and ileum, and the amount of drug in the gut lumen (A i ) is controlled by the gut transit rate constant (K t,i ) and absorption rate constant (K a,i ). The change in mass within each gut lumen was described as follows (Qian et al., 2019): where i 1, 2, and 3 corresponds to the duodenum, jejunum, and ileum, respectively. The calculation of the K a,i value is given by the following: 4 where K a,i0 represents the initial absorption rate constant and f a represents the absorption adjustment factor. The calculation of the K a,i0 value is given by the following: 6 K a,i0 2 × P eff,human × ASF r i , where P eff, human is the effective permeability in humans and r i is the average radius of different intestinal segments (the duodenum, jejunum, and ileum). ASF is the absorption scale factor.

Enterocyte Compartment
Each different intestinal segment is linked with its enterocytes compartment. The amount of drug in each enterocyte compartment (A ent,i ) was calculated according to the following equations (Qian et al., 2019): where Q ent,i and V ent,i represent the blood flow rate and the volume of each enterocyte compartment, respectively. V pv is the volume of the portal vein compartment, f ugut is the free drug concentration in the enterocyte compartment, and K p,ent is the ratio of intestine over blood drug concentration. CL int,ent (t) and Abundance ent,i (t) represent the intrinsic metabolic clearance and CYP3A content in the enterocyte compartment with inhibitors/inducers or without inhibitors/ inducers, respectively. V max and K m are metabolic parameters of the drug in each different enterocyte compartment. According to the reported result, saxagliptin is a weak P-glycoprotein substrate (Boulton, 2017); consequently, the efflux effect of P-glycoprotein on the amount of saxaliptin in the gut lumen was not considered in this study.

Portal Vein Compartment
The drug from the enterocyte compartment enters the liver via the portal vein, and hence, the change of drug amount over time in the portal vein compartment (A pv ) is illustrated according to the following equations (Li et al., 2012;Qian et al., 2019): The PBPK-DO model includes the gut lumen, enterocytes, the portal vein, the liver, blood (arterial supply and venous return), the lung, and other eliminating tissues (kidney) and non-eliminating tissues (11 compartments). The CYP3A metabolism enzyme was set into enterocytes and liver compartments. DPP-4 was assumed to only reside in the venous return compartment. 4 PBPK book.pdf 6 GastroPlusManual-9.7.pdf Frontiers in Pharmacology | www.frontiersin.org October 2021 | Volume 12 | Article 746594 where C ab is the drug concentration in arterial blood.

Liver Compartment
The liver is a main eliminating tissue of saxagliptin, and the amount of drug in the liver is described as follows 6 (Li et al., 2012): where A li , V li , Q li , and Q lia are the drug amount in the liver, the volume and hepatic blood flow, and the hepatic artery blood flow rate of the liver, respectively. Rbp is the blood-to-plasma concentration ratio, f up is the fraction of free drug in the plasma, and K p,li is liver-to-plasma partition coefficient. M is molecular weight of the drug. CL int,li (t) represents the intrinsic metabolic clearance of a drug in the liver. Abundance li (t) represents the hepatic CYP3A amount with inhibitors/inducers or without inhibitors/ inducers. ISEF is the intersystem extrapolation factor. Owing to the finding that saxagliptin is not almost eliminated through biliary excretion (Fura et al., 2009), biliary clearance of saxagliptin is not incorporated into this PBPK model.
The amount of metabolite in the enterocyte compartment (A met,ent ) is illustrated by the following: The amount of metabolite in the liver (A met,li ) is illustrated by the following: where f scale is the scale factor of metabolite conversion. The total amount of metabolite (A met ) enters the venous blood compartment directly without consideration of gut absorption and is calculated by the following: where K a,met0 and K a,met are the initial absorption rate constant and absorption rate constant of metabolite M2, respectively. f a,met represents the absorption adjustment factor of the metabolite.

Kidney Compartment
The change of drug amount with time within the kidney (A ki ) is described by the following (Li et al., 2012): where Q ki and V ki are the blood flow and the volume of the kidney, respectively. K p,ki is the kidney-to-plasma partition coefficient. CLr is the kidney elimination rate of the drug.

Lung Compartment
The amount of drug in the lung (A lu ) is described as follows (Li et al., 2012): where Q lu and V lu are the blood flow and the volume of the lung, respectively. K p,lu is lung-to-plasma partition coefficient.

Other Non-eliminating Tissue Compartments
Overall, the amount of drug within non-eliminating tissue compartments (A nt ) (the adipose, the bone, the brain, the heart, the muscle, the skin, the spleen, the red marrow, the yellow marrow, reproductive organs, and the rest of the body), except for the lung, follows the following equation (Li et al., 2012): where Q nt and V nt are the blood flow sum and the mean volume of other non-eliminating tissues, respectively. K p,nt is the partition coefficient sum of other non-eliminating over plasma.

Arterial Blood Compartment
The drug concentration within the arterial blood (C ab ) is described as follows (Li et al., 2012;Fu et al., 2019): where V ab is the volume of the arterial blood compartment. i represents all eliminating and non-eliminating tissues except for the lung tissue.

Venous Blood Compartment
The drug concentration within the venous blood (C vb ) is described as follows (Li et al., 2012;Fu et al., 2019): where V vb is the volume of the venous blood compartment. i represents all eliminating and non-eliminating tissues except for the lung tissue.

CYP3A Dynamics of Inhibition and Induction
Here, the inhibition and induction model of CYP3A were developed assuming that intestinal inhibition and induction dynamic parameters were identical to hepatic dose. The CYP3A reversible competition inhibition dynamics is described as follows (Baneyx et al., 2014): The CYP3A time-dependent inhibition dynamics is described as follows (Baneyx et al., 2014): The CYP3A induction dynamics is described by the following (Qian et al., 2019): where I li is the concentration of the inhibitor/inducer within the targeted tissue (the gut and liver). f up,in is the fraction of the free inhibitor/inducer in the plasma. K p,in is the partition coefficient of the inhibitor/inducer of the targeted tissue (the gut and liver) to plasma. M in is the molecular weight of the inhibitor/inducer. K i and k inact are inactivation parameters of the inhibitor against CYP3A4. EC max and EC50 are inductive parameters of the inducer on CYP3A4. k deg is the degradation rate constant of CYP3A4, and in this study, it was assumed that the k deg value in the liver is identical to that in the gut.

DPP-4 Engagement Dynamics by Saxagliptin and Metabolite M2
The DPP-4 engagement time course by saxagliptin and metabolite M2 is described by a series of equations as follows (de Witte et al., 2016): where i 1 and 2 corresponds to saxagliptin and active metabolite M2, respectively. TC is the concentration of saxagliptin/M2-DPP-4 complex formed. T free is the concentration of free DPP-4. T total is the sum of TC plus T free . C vb,i , k on,i , and k off,i are the drug concentration within the venous blood and the on-rate and off-rate of saxagliptin and active metabolite M2, respectively.

DDI Prediction Through the PBPK-DO Model
The DDI predictions were conducted to evaluate the parameter changes of PD (TO AUC , area under the occupancy-time curve; TO max , maximum occupancy; and DTO >60% duration of >60% TO) and PK (AUC, area under the concentration-time curve; and C max , peak concentration) on saxagliptin and M2 with and without inhibitors/inducers.

Development of the PBPK-DO Model for Saxagliptin and Metabolite M2
The PBPK-DO model for saxagliptin has been established by a massive number of input parameters from Table 1 and Supplementary Table S4. Figures 2A,B show the predictions and observations of the human PK and DPP-4 occupancy for saxagliptin and metabolite M2, respectively, after oral administration of a 5 mg saxagliptin. Comparison of observed PK and DPP-4 occupancy parameters with simulated parameters of saxagliptin and M2 is summarized in Supplementary Table S5. It is clearly indicated that human PK simulation corresponds closely to observed values for saxagliptin and M2 (Upreti et al., 2011) and that simulation of time course of DPP-4 occupancy in humans by saxagliptin could also be matched with experimentally determined values very well (Upreti et al., 2011). The simulation results have displayed that the developed PBPK-DO model could accurately predict PK profiles and DPP-4 time profiles in humans for saxagliptin and M2.

Effect of Inhibitors/Inducers on PK and PD of Saxagliptin in Humans
The PK profile predictions of inhibitors/inducers have been shown in Supplementary Figure S1, and the predicted and observed data are listed in Supplementary Table S6. The accuracies of PK prediction using the developed PBPK model have been verified by comparing predicted and observed PK of the inhibitors/inducers. The comparison displayed that predicted PK profiles matched observed profiles well (Supplementary Figure S1) and that all fold errors were less than 2 between predicted and observed PK data (Supplementary Table S6). The result indicated that the established PBPK model could accurately simulate the PK process in humans of inhibitors and inducers.
The PK profiles and DPP-4 occupancy profiles of saxagliptin in humans at a dose of 5 mg were simulated using the developed PBPK-DO model after coadministration of ketoconazole (200°mg, twice daily), delavirdine (400°mg, twice daily), and rifampicin (600°mg, once daily) for 10°days, respectively (Figures 3, 4). The C max and AUC 0-t of saxagliptin coadministrated with multidose ketoconazole on the sixth day increased by 1.89-and 3.42-fold of those without ketoconazole, respectively, (Supplementary Table S7). Slight underestimation was found in C max and AUC of saxagliptin for DDI predictions (C max ratio: 1.89; AUC 0-t ratio: 3.42) versus DDI observations (C max ratio: 2.44; AUC 0-t ratio: 3.67) after coadministration of ketoconazole 7 . It was observed that C max and AUC of saxagliptin increased by 1.33-and 1.50-fold with delavirdine compared to those without the inhibitor, respectively (Supplementary Table S8). The PK profiles and parameters of saxagliptin were strongly influenced following coadministration of rifampicin, while PK profiles and parameters of metabolite M2 were nearly unchanged ( Figure 3D and Supplementary Table S9). Slight overestimation was observed in C max and AUC of saxagliptin compared to actually experimentally determined values, with DDI predictions of a C max ratio of 0.31 and an AUC 0-t ratio of 0.19 versus DDI observations of a C max ratio of 0.42 and an AUC 0-t ratio of 0.24 (Upreti et al., 2011). The predicted AUC 0-t ratio of M2 was slightly below clinical experiment data (0.78 versus 0.91) (Upreti et al., 2011), while the C max ratio of M2 had medium differences between predicted and experimentally determined values (0.88 versus 1.38) (Upreti et al., 2011).
The levels of DPP-4 occupancy by saxagliptin (5 mg) and M2 coadministrated with multidose ketoconazole have been significantly enhanced compared to that without ketoconazole, with the lowest DPP-4 occupancy being >80% ( Figure 4A, red line). Owing to significant improvement of DPP-4 occupancy, next, the time course of DPP-4 occupancy was simulated at a lower dose of 2.5 mg ( Figure 4A, green line). The comparison between the DPP-4 occupancy profile at 5 mg of saxagliptin without ketoconazole (black line in Figure 4A) and that at 2.5 mg of saxagliptin with ketoconazole (green line in Figure 4A) demonstrated that oral saxagliptin should be decreased to 2.5 mg when there is coadministration with ketoconazole. The result was in good agreement with reported data in the literature 7 . Although AUC of saxagliptin increased 1.5-fold with delavirdine, however, the DPP-4 occupancy time profile by saxagliptin and M2 coadministrated with multidose delavirdine almost coincided with that without delavirdine ( Figure 4B). In contrast, while C max and AUC of saxagliptin with coadministration of rifampicin had a considerable decrease, the percent occupancy of DPP-4 by saxagliptin and M2 with coadministration of rifampicin closely resembled that without rifampicin ( Figure 4C). This prediction was in line with the published result in the study (Upreti et al., 2011). The similar occupancy of DPP-4 indicated that M2 was likely the main contributor to human DPP-4 occupancy under the circumstance of coadministration with rifampicin.

DISCUSSION
In this study, we have developed a PBPK-DO mathematical model characterized by two simultaneous simulations (parent/metabolite and PK/PD), which was first utilized to quantify the impacts of CYP3A4 FIGURE 3 | Human plasma level of saxagliptin and M2 following coadministration of inhibitors or inducers. Human plasma profiles of saxagliptin with ketoconazole (A, 200°mg, twice daily), with delavirdine (B, 400°mg, twice daily) and rifampicin (C, 600°mg, once daily) and of M2 with rifampicin (D). The blue squares, blue up-triangles, and red uptriangles refer to observed PK data of saxagliptin without inhibitors/inducers (□) and observed PK data of saxagliptin (△) and M2 (△) with rifampicin on the sixth day, respectively. inhibitors and inducers on the PK and PD of saxagliptin and M2 in humans simultaneously. Ketoconazole and rifampicin are recommended as the standard CYP3A4 competitive inhibitor and inducer for the potential clinical DDI study, respectively. In addition, we also assessed the effect of a TDI (delavirdine) on the PK and PD of saxagliptin and M2, which is a distinct inhibition type from ketoconazole. Although some studies have showed that P450 metabolic activity of a drug could be frequently different in the gut and the liver (Choe et al., 2017;Kapetas et al., 1208;Kurucz et al., 2019;Klomp et al., 2020), the metabolic parameter of saxagliptin by intestinal CYP3A4 was considered comparable to hepatic CYP3A4 in this PBPK-DO model. f a and ASF were incorporated into this model to optimize PK peak time of saxagliptin and scale the effective permeability 4, 6 , which greatly improved prediction performance. According to the literature (Upreti et al., 2011;Boulton, 2017), here, we assumed that all DPP-4 enzymes were located in the blood compartment in the present model. Across published studies (Rowland Yeo et al., 2011;Bolleddula et al., 2021), k deg values of CYP3A4 ranged from 0.0077 -1 to 0.03 h −1 . However, a recent study has confirmed that it was a more reasonable value at 0.03 h −1 for the most accurate prediction (Rowland Yeo et al., 2011); hence, the k deg value was set at 0.03 h −1 in this PBPK-DO model. The induction activity (E max and EC 50 ) of rifampicin suggested major individual variability between different literature studies (Yamazaki et al., 2015;Almond et al., 2016;Asaumi et al., 2018;Hanke et al., 2018). To minimize the variation, we used mean values of E max (6.2 μm) and EC 50 (0.6 μm) in the current model (Qian et al., 2019). Use of rifampicin induction parameters could have been rationalized indirectly that the predicted expression amount of mean intestinal CYP3A4 after induction by rifampicin was close to the experimentally determined value in humans 8 (6.4-fold change versus 4.4-fold change). The content of human microsomal protein was set at 38 mg/g in the liver 6 rather than 45 mg/g liver in some studies (Guo et al., 2013;Qian et al., 2019) because the former parameter had been proved to be a more reasonable value in more studies and built-in business software 6 .
All predictions were within 2-fold of observations among them, and the highest fold error was 1.5-fold, which occurred between the predicted and observed C max ratios of M2 with rifampicin. In accordance with FDA clinical DDI guidance (Sudsakorn et al., 2020), if the AUC ratio of a drug with and without the inhibitor is ≥1.25 or ≤0.8 with and without the inducer, clinically relevant DDI should be considered. In our simulations, the AUC ratio of saxagliptin was found to be 1.50 with delavirdine and to be 0.19 with rifampicin, which could occur in significant clinical DDIs based predominantly on PK comparisons. Nevertheless, DPP-4 occupancy by saxagliptin with delavirdine or rifampicin was almost unchanged, and PD simulations displayed that delavirdine or rifampicin would not cause human DDIs for saxagliptin, which was in good agreement in clinical experiments (Upreti et al., 2011). The simulation results also have further demonstrated the importance of two simultaneous simulations in the present model.
It was reported that the C max and AUC 0-t of saxagliptin increased by less than 2-fold in patients with severe hepatic impairment but by more than 2-fold in patients with severe renal impairment, respectively (Boulton, 2017). PK simulation in humans with hepatic or renal impairment is performed by modulating many physiological parameters based on healthy humans (Malik et al., 2020). However, currently, this simulation model cannot simulate PK of patients with hepatic or renal impairment yet. Recent literature reported that catalytic activities of 27 CYP3A4 variants on the in vitro metabolism of saxagliptin were evaluated (Liu et al., 2021). CYP3A4 variants showed decreased activities ranging from 1.9 to 77.1% as compared to the wild type. Hence, we also preliminarily evaluated the effect of genetic variations in metabolizing enzymes on PK and PD of saxagliptin using this model. Here, we only simulated the PK and PD of saxagliptin in humans with variant CYP3A4*22. In the simulation, expression and activity of CYP3A4 were replaced with 59 and 40% of the wild type (Alqahtani and Kaddoumi, 2016). The results displayed that C max and AUC 0-t of saxagliptin increased by about 2-fold, and DPP-4 occupancies by saxagliptin between three oral doses (5, 2.5, and 1 mg) have a slight difference (Supplementary Table S10 and Supplementary Figure S2). Of note, due to clinical data unavailability, prediction accuracy need be proven with further in vivo studies.

CONCLUSION
Taken together, this mathematic model is characterized by two simultaneous simulations (parent/metabolite and PK/PD), describing two interaction processes between inhibitors/inducer-CYP3A4 and saxagliptin/M2-DPP-4. We conceive that compared to most current single PK-DDI predictions, the wide application of the PBPK-DO model has the power to improve the predictions of potential clinical DDIs for victim drugs metabolized by CYP3A4.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
Conceptualization, HL, YL, and GW; methodology, GL, BY, and GW; software, GL and FP; validation, BY and JL; formal analysis, GW, GL, and BY; investigation, JL; resources, XJ; data curation, WY; writing-original draft preparation, GW, GL, and BY; writing-review and editing, GW, GL, and BY; visualization, GL and BY; supervision, GL and BY; project administration, HL, YL, and GW. All authors have read and agreed to the published version of the article.